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Abstract 

We investigate the dynamics of relaxation, by surface tension, of a family of curved interfaces between an inviscid and 
viscous fluids in a Hele-Shaw cell. At t = the interface is assumed to be of the form \y\ = A x m , where A > 0, m > 0, 
and x > 0. The case of < m < 1 corresponds to a smooth shape, m > 1 corresponds to a cusp, whereas m = 1 
corresponds to a wedge. The inviscid fluid tip retreats in the process of relaxation, forming a lobe which size grows 
with time. Combining analytical and numerical methods we find that, for any m, the relaxation dynamics exhibit 
self-similar behavior. For m 7^ 1 this behavior arises as an intermediate asymptotics: at late times for < m < 1, 
and at early times for m > 1. In both cases the retreat distance and the lobe size exhibit power law behaviors in 
time with different dynamic exponents, uniquely determined by the value of m. In the special case of m — 1 (the 
wedge) the similarity is exact and holds for the whole interface at all times t > 0, while the two dynamic exponents 
merge to become 1/3. Surprisingly, when m^l, the interface shape, rescaled to the local maximum elevation of the 
interface, turns out to be universal (that is, independent of m) in the similarity region. Even more remarkably, the 
same rescaled interface shape emerges in the case of m = 1 in the limit of zero wedge angle. 
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1. Introduction 

Consider a curved interface between a low- 
viscosity fluid (for example, water) and a high- 
viscosity fluid (for example, oil) in a large horizon- 
tal Hele-Shaw cell. If the system is unforced, the 
interface undergoes relaxation by surface tension, 
ultimately approaching either a straight line, or 
a circle (or breaking into several domains which 
then become circles). This free boundary problem 
is hard to solve, as the governing equations (see be- 
low) are non-local. This is especially true when the 
interface initially has a complex shape, as observed 
in numerous experiments when the viscous fluid is 
displaced by the inviscid fluid in radial geometry, 
see Ref. [1] and references therein. The initial shape 



complexity is caused here by the viscous finger- 
ing instability that develops during the preceding 
forced stage of the flow [2,3]. The (strongly) forced 
Hele-Shaw flow is a standard paradigm in fluid 
dynamics and nonlinear dynamics [4,5,6,7]. The 
(small) surface tension there is usually invoked in 
order to regularize the otherwise singular dynamics 
on small scales. We are interested in this paper in 
an unforced flow, where surface tension is the only 
driving mechanism. Here is the formulation of the 
unforced Hele-Shaw (UHS) flow model that we will 
be dealing with throughout this paper. Let one 
fluid have negligible viscosity, so that the pressure 
in it is uniform. The velocity of the viscous fluid is 
v (r, t) = -(H 2 /I2fi) Vp (r, t), where p is the pres- 
sure, fj, is the dynamic viscosity, and H is the plate 
spacing [2,3,4,5]. Therefore, the interface speed is 
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v n = -(H 2 /12^)d nP , (1) 

where index n denotes the components of the vec- 
tors normal to the interface and directed from the 
inviscid fluid to the viscous fluid, and d n p is evalu- 
ated at the corresponding point of the interface 7. 
The viscous fluid is incompressible. Therefore, its 
pressure is a harmonic function: 

V 2 p = 0. (2) 

The Gibbs-Thomson relation yields a boundary con- 
dition at the interface: 

p | 7 = (tt/4) oK. , (3) 

where a is surface tension, and IC is the local cur- 
vature of the interface, positive when the inviscid 
region is convex outwards. Finally, as the flow is un- 
forced, we demand 

Vp=0 at r^oo. (4) 

Equations (l)-(4) define the UHS problem (see Rcfs. 
[8,9,10] and references therein for a more detailed 
discussion). The UHS model gives an instructive ex- 
ample of non-local area-preserving curve-shortening 
dynamics. 

The UHS flow (l)-(4) is not intcgrable. Moreover, 
until recently even no particular analytic solutions 
to this class of flows had been found, except for the 
simple solutions provided by a linear stability anal- 
ysis of a single, slightly deformed flat or circular in- 
terface [11]. Recently, some analytic solutions have 
been obtained for two special initial interface shapes. 
In the first of them, the inviscid fluid domain at t — 
has the form of a half-infinite stripe [12]. As time 
progresses, the tip of the stripe retreats and develops 
a lobe. At long times, the growing lobe approaches a 
self-similar shape, whereas the lobe size and retreat 
distance follow a power law in time with different 
dynamic exponents: 1/5 and 3/5, respectively. 

In the second case the assumed form of the invis- 
cid fluid was a wedge [13]. As this initial condition, 
and Eqs. (l)-(4), do not introduce any length scale 
into the problem, the solution is self-similar at all 
times, with a single dynamic exponent 1/3 [13]. The 
scale-invariant interface shape in this case is given 
by the solution of an unusual inverse problem of po- 
tential theory. Gat et al. [13] solved this problem 
perturbatively for an almost flat wedge, and numer- 
ically for several values of the wedge angle. 

The results of Refs. [12,13] suggest that the values 
of dynamic exponents, and other attributes of the 



self-similar asymptotics, are determined by the ini- 
tial shape of the retreating edge of the inviscid fluid 
domain, while the two solutions obtained in Refs. 
[12,13] are particular members of a broader family 
of solutions. The results of the present paper con- 
firm this scenario. We consider here a more general, 
power-law shape \y\ — Ax m , where A > 0, m > 0, 
and x > 0, and show that, for any m > 0, the re- 
laxation dynamics exhibit self-similar intermediate 
asymptotics: a late-time asymptotics for < m < 1 
and an early-time asymptotics for m > 1. The re- 
treat distance and the lobe size show, for any m, 
a power law behavior in time with exponents and 
pre-factors uniquely determined by m. The case of 
m = 1, investigated in Rcf. [13], is special: here the 
self-similarity is exact and occurs for all times t > 

0. while the two dynamic exponents merge and be- 
come equal to 1/3. Surprisingly, at m =/= 1 the in- 
terface shape, rescaled to the local maximum eleva- 
tion of the interface, turns out to be universal (that 
is, independent of m) in the similarity region. Even 
more remarkably, the same rescaled interface shape 
emerges in the case of m = 1 in the limit of zero 
wedge angle. 

Here is a layout of the rest of the paper. In Section 
II we generalize to an arbitrary m > the approach, 
suggested by Vilenkin et al. [12] for a half-infinite 
stripe (m = 0). We present there a simple asymp- 
totic scaling analysis that predicts (i) the dynamic 
exponents of the self-similar part of the flow, (ii) the 
exponent of the power-law tail of the scale-invariant 
shape function of the interface, and (iii) the validity 
range of the scaling behavior at < m < 1 and m > 

1. In Section III we report the results of a numerical 
solution of (l)-(4) for the cases of m = 1/4, 1/2 and 
5/4, compare them with our theoretical predictions 
and report the universality of the rescaled interface 
shape. In addition, we present in Section III our new 
numerical results for small-angle wedges. A brief dis- 
cussion and summary are presented in Section IV. 

2. Interface dynamics: theoretical 
predictions 

Let at t = the interface shape be \y\ = Ax m , 
where A > 0, m > 0, and x > 0. The case of < 
m < 1 corresponds to smooth shapes, m > 1 cor- 
responds to a cusp, see Fig. 1, while m = 1 cor- 
responds to a wedge. The parameter A has the di- 
mension of length 1 "™ 1 . This implies that the case of 
m = 1, investigated by Gat et al. [13], is special, 



2 



a) 

f 


r 


L 


— "^"^ 

/ y 
/ y 
/ y 
1 *»» 

'"" X 


b) 


\ "Si 

\ V. 

\ V. 



Fig. 1. A schematic setting for the interface dynamics for 
< m < 1 (a) and m > 1 (b). 

as the parameter A does not introduce any length 
scale there. We assume in this Section that to ^ 1 
and measure the distance in units of A = yl 1 /( 1_m ) ; 
the time in units of t = 48 /i A 3 /(nail 2 ), and the 
pressure in units of po = 7rcr/(4A). In the rescaled 
variables the interface shape is |y| = x m , while Eqs. 
(1) and (3) are parameter- free: 



and 
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~d n p 



JC. 



(5) 



(6) 



Because of the finite surface tension, the pres- 
sure gradient in the viscous fluid is largest at the 
tip, so the tip will retreat along the x-axis and de- 
velop a lobe which size will grow with time. Let us 
assume (and then test numerically) that the evolv- 
ing interface can be fully characterized by two time- 
dependent length scales: the curvature radius, R(t), 
and the retreat distance, L(t), of the tip. In the lead- 
ing order this assumption implies a similarity ansatz, 
in the moving frame x\ = x — L(t), for the interface 
elevation y — h(x, t) in the lobe region and region II 
(see below): 



h s (x u t) = R(t)<S> 



X\ 



R(t) 



(7) 



In the scaling regime, the dynamic length scales L(t) 
and R(t) exhibit a power law behavior: 

L(t)=at a and R(t) = bt f) , (8) 

where the exponents a and (3 and coefficients a and 
b are m-dependent. The objectives of this work is 
to determine analytically a, [3 the exponent of the 
power-law tail of the scale-invariant shape function 
of the interface <& \x\/R(t)\, and the validity range 
of the scaling behavior at < m < 1 and to > 1. 
These analytic findings will be verified and comple- 
mented by our numerical solutions which, in addi- 
tion, provide the whole scale-invariant shape func- 
tion <!> [x\/R(t)] and the coefficients a and b of the 
power laws of R(t) and L(t) for four different values 
of m, in Sec. III. 

A simple asymptotic scaling analysis for m ^ 1 
is possible because of the following property of the 
initial shape of the inviscid fluid region: it has a flat 
region at cither very large distances i»l (for < 
m < 1), or very small distances x <g; 1 (for m > 1). 
We will see that, as a result, the scaling regime will 
hold at t > 1 for < m < 1 and at t < 1 for 
m > 1. The calculations are very similar to those 
of Ref. [12], where the particular case of m = 
was considered. Introduce for a moment the system 
of polar coordinates r and 4> with the origin at the 
moving tip of the interface, as shown in Fig. la. In 
view of the Gibbs-Thomson condition (6), p must 
vanish, in the leading order, at the flat region, that 
is at and <fi —> 27r. This corresponds either to 
the region 1 <C R(t) <C r for < m < 1, or to the 
region R(t) <C r <C 1 for to > 1. Furthermore, in 
view of the same Eq. (6), p = K, ~ 1/R(t) in the lobe 
region (for definiteness, at (j> = ±7r/2). Therefore, 
the leading term of the far-field multipole expansion 
of p [14] can be written as 

p(r,ct>,t) = C m [R(t)r}- 1/2 sin(0/2), (9) 

where C m = 0(1). Having demanded the boundary 
condition (6) in Eq. (9), we stretched the validity 
region of Eq. (9), r » R(t), toward r ~ R(t), but 
this can only affect the value of constant C m that 
our theory cannot give anyway. Now we can employ 
Eqs. (5) and (9) and estimate the normal component 
of the interface speed in the far field region. For the 
upper interface of the far field region 

Vn = -\% { ^ Q) = - 2R^{t)rm - ^ 

At this point we return to the Cartesian coordinates 
xi,y. In the far field region, that is at 1 <C R(t) <C 
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x\ for < to < 1 or R(t) <C x\ -C 1 for m > 1, the 
quantity dh{x\,t)/dt is given by Eq. (10), and we 
obtain 



r 

-x m = J 



Cm dt' 



2R l / 2 {t')x[ 



3/2 



(11) 



As we will check a posteriori, L(t) is always much 
greater than R{t) in the scaling regime. How can we 
simplify the calculation of the intergal in Eq. (11)? 
The far field region x\ R(t) can be divided into 
two sub-regions: x\ » L{t) (region I) and R(t) <C 
xi <C L{t) (region II, or neck). In region I we can, in 
the leading order, take x\^ 2 out of the integral and 
arrive at 



h(xi,t) 



x ~ 



o 3/2 

2x{ 



dt' 



t 



3/2 



i?V2(t) 



,(12) 



where we have assumed that R{t) is a power of t, 
and disregarded the coefficient C m = 0(1). 

In region II at fixed x, the main contribution to 
the integral in Eq. (11) comes from times close to 
t, so that xi(t)/L(t) <^Lt — t' < t. Indeed, we can 
expand xi(t') = xi(t) + L(t)(t — t') + . . . and, in 
the leading order, neglect higher order terms. The 
effective time interval for the integration is (t—St 1 , t), 
where St' ~ x\ (t)/L(t). Furthermore, R}l 2 {t') can 
be evaluated at t' — t, as its variation on the time 
interval (t — St', t) is negligible. Then, extending the 
lower limit of the integral to — oo and calculating the 
remaining elementary integral, we obtain 



As, by assumption, L(t) is a power law, SAi is com- 
parable to 5 An. As we will check shortly, the con- 
tribution to the area of the lobe region itself, SAr ~ 
R 2 (t), is negligible compared to SAi and 8 An as 
long as we are in the scaling regime (t ^> 1 for < 
to < 1 or t <C 1 for m > 1). 

Now we employ the exact integral of motion of the 
system: the area conservation of each of the fluids. 
The area loss because of the retreat [which is equal 
to L{t) m+1 /{m + 1) - L{t) m+1 ] must be equal, in 
the leading order, to the area gain in the far field 
region. This follows 

L(t) m+1 ~ SAi(t) ~8A n {t), 

which yields a relation between the two dynamic 
length scales R(t) and L(t). Another relation be- 
tween these two quantities follows from Eq. (9) . We 

obtain Vi dp/dr [r ~ R(t), (f> ~ n] - R~ 2 (t), 

and demand 

L(t) - R~ 2 {t) . 

These two relations immediately yield the dynamic 
exponents a and [3: 



a = 



h(xi,t) 



C m 



1/2 ' 



(13) 



where the factor C m — 0(1) is in excess of accuracy 
and can be disregarded. Now we can estimate the 
contributions of regions I and II to the area gain SA 
in the far field region. In region I (correspondingly 
II) the main contribution to the integral over x\ 
comes from the lower (correspondingly, upper) limit 
of integration. Therefore, we integrate Eq. (12) over 
x\ from, say 2L(t) to infinity and Eq. (13) from 
R{t) to 2L(t). The results are: 



SAi(t) ^ 

and 

SAnit) 



t 



LV2( t ) R}/ 2 {t) 
^ L{t)R^I\t) 



in region I , 



in region II . 



(14) 



(15) 



4to + 5 



and (3 



2m + 1 
4to + 5 



(16) 



Once the scaling relations for L(t) and R(t) are 
found, we can calculate [up to m-dependent numer- 
ical pre-factors 0(1)], additional quantities. For ex- 
ample, the interface elevation in region I becomes 

h(x, t) - x m - ti|Stj x - 3 / 2 , 
see Eq. (12). In region II, see Eq. (13), we obtain 



h(x, t) 



3(2m + l) 
£2(4 m + 5) x - 



1/2 



Importantly region II belongs to the similarity re- 
gion, as was first observed in Ref. [12] in the case of 
to = 0. In this region <&(£) ~ £ -1 / 2 : a universal (m- 
independent) power law of the similarity variable 
£ = xi/R(t). The presence of the decreasing asymp- 
tote $(£) ~ £ -1 / 2 implies that, for any m ^ 1, the 
shape function $(£) must have a local maximum. 

Once we obtained the solution, we can check it for 
self-consistency with all the assumptions we made. 
First, it can be easily checked that, in the scaling 
regime t ^> 1 (for < m < 1) or f < 1 (for 
to > 1) we have L(t) ^> R(t), as we assumed. Now, 
the lobe area SAr ~ R 2 (t) grows with time as ~ 
t 2(2m+i)/(4m+5) Tnig value is j^eed much less than 

<L4/(i) - SA H (t) - f3(m+l)/(4m+5) j n thc sca l ing 

regime. Furthermore, thc time-dependent interface 
elevation in the lobe region R(t) ~ L(i)( 2m+1 )/ 3 is 
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much larger, in the scaling regime, than the initial 
interface elevation at the moving tip, which is L(t) m . 
This strong inequality serves as a necessary condi- 
tion for the assumed similarity asymptotics (7) as a 
leading order description. 

3. Numerical solution 

In order to test the predicted dynamic exponents 
a and (3 and verify the presence of the self-similar 
regime, we solved the problem numerically for three 
different values of parameter to. In addition, we re- 
turned to the case of to = 1, previously considered 
in Ref. [13], and solved it numerically for small val- 
ues of the wedge angle. 

3.1. Numerical method 

Our numerical algorithm [16] employs a variant of 
the boundary integral method for an exterior Dirich- 
let problem formulated for a singly connected do- 
main, and explicit tracking of the contour nodes. 
The harmonic potential is represented as a potential 
produced by a dipole distribution with an a priori 
unknown density D on the contour. The dipole den- 
sity D is found numerically from an integral equa- 
tion which is a modification of the well-known jump 
relation of the potential theory [17]. Computing an- 
other integral of this (already found) dipole density, 
one obtains the harmonic conjugate function, whose 
derivative along the contour yields, by virtue of the 
Cauchy theorem, the normal velocity of the inter- 
face. 

We used a piecewise constant function for a dis- 
crete approximation of D and a piecewise linear 
function for discretizing the interface. The integral 
entering the integral equation is represented as a 
sum of D multiplied by a kernel which is integrated 
analytically between two neighboring nodes. This 
approximation was previously suggested for the in- 
ner problem [15]. We found that it is also efficient 
in the outer problem in the following cases: (i) for 
a long and slender domain, (ii) in the vicinity of 
the cusp at to > 1, and (iii) for a wedge (to = 1) 
with a small angle. The numerical approximation 
is described in detail in Ref. [16]. The method re- 
quires an inhomogeneous grid with a small spacing 
in regions of high curvature of the contour, and we 
used a grid with spacing exponentially growing with 
the distance from the interface's tip. The number of 
grid nodes was reduced as the interface's perimeter 



decreased, and the curvature radius of the tip in- 
creased. 

The shape of the numerical interface at t = is de- 
termined by the following parameters: the exponent 
to, the domain size A > 0, and the cutoff parameter 
e > that was used for to > 1, see below. One quar- 
ter of the interface is represented as a graph h(x) = 
(x + A) m , where —A + e < x < 0. The second quar- 
ter is obtained by reflecting this graph with respect 
to the x-axis. Then, by reflecting the two branches 
with respect to the y-axis, we obtained the closed 
interface we worked with. In this manner we could 
exploit the four-fold symmetry of the domains and 
achieve a four-fold reduction in the number of al- 
gebraic equations, approximating the integral equa- 
tion. For to < 1, the tip is smooth, and we took e = 
0. For to > 1 there is a cusp at x = — A that our nu- 
merical method can not handle. A similar difficulty 
arises for m = 1 if the wedge angle is very small. A 
positive e allows one to employ the method, if the 
node spacing in the cutoff region is less than e m . 

We measured the retreat distance of the tip L(t) 
and the interface shape at different rimes for three 
values of m: m = 1/4, 1/2 and 5/4, and also for 
three different wedge angles for m = 1. In the pro- 
cess of relaxation, each of these domains ultimately 
becomes a perfect circle. Therefore, to observe the 
self-similar asymptotics we performed the measure- 
ments at times much shorter than the characteristic 
time of relaxation toward a circle. In addition, to 
minimize the influence of other tips, we performed 
the measurements sufficiently close to a chosen tip 
(at distances much smaller than the distance be- 
tween the chosen tip and the neighboring tip) . These 
two limitations are especially relevant at to < 1, 
where theory predicts self-similarity at sufficiently 
long times. On the contrary, at m > 1 we performed 
the measurements at very short times. Here the main 
limitation comes from the presence of the cutoff, 
which necessitates a sufficiently long "waiting time" 
so that the influence of the cutoff on the solution 
can be neglected. 

For to = 1/4 we took A = 2 x 10 4 and the mini- 
mum node spacing SS = 1/2. The spacing increased 
exponentially with the distance from the tip, and 
the initial number of nodes was 2 x 10 3 (here and 
in the following - per quarter of the interface). For 
to = 1/2 we chose A = 10 5 , SS = 1 and the initial 
number of nodes 10 3 . The set of numerical param- 
eters for to = 5/4 was A = 1CT 3 , SS = 1(T 8 , e = 
(2 x 1(T 8 ) 4/5 = 6.93 x 1(T 7 , and the initial number 
of nodes 1200. 
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For m = 1 and the wedge angles 9 = 10° and 5° 
we used A = 1, SS — 10~ 5 , and the initial number 
of nodes 10 3 . The cutoff parameter e was chosen to 
be e = 3SS cot 6 « 0.017 for 6 = 10°, and 0.034 for 
6 = 5°. For = 2° we took A = 3 • 10" 3 , SS = lO" 6 , 
the initial number of nodes 2 ■ 10 3 , and the cutoff 
parameter e w 2.8 • 10 -5 . 

In all cases the time step was chosen to be 10 ~ 3 
times the maximum of the ratio of the interface cur- 
vature radius and the interface speed at the same 
node. This choice of numerical parameters was dic- 
tated by the fact that, at < m < 1, we were in- 
terested in a long-time behavior, whereas at m > 1 
we needed to focus on very earlier times, in order to 
observe the predicted self-similarity and scalings. 
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Fig. 2. Snapshots of a part of the simulated system for 
m = 1/4 at times t = (the black solid line), 2.7 X 10 5 (the 
red dashed line) and 1.2 X 10 6 (the blue dotted line). No- 
tice the large difference between the horizontal and vertical 
scales. 

3.2. Numerical results 

Figure 2 shows snapshots of a part of the sys- 
tem for m = 1/4 at times t = 0, 2.7 x 10 5 and 
1.2 x 10 6 . One can see a lobe developing and growing 
with time. Shown in Fig. 3a is the retreat distance 
L{t) versus time. A power law fit yields exponent 
0.50 which coincides with the predicted theoretical 
value 1/2, see Eq. (16). It is more convenient numer- 
ically to measure the local maximum height of the 
interface h max (t), rather than the curvature radius 
at the tip R(t). Because of the self-similarity, the 
quantities h max (t) and R(t) are expected to exhibit 
the same power law dependence (of course, with dif- 
ferent pre- factors). Fig. 3b shows h max versus time 
in the case of m = 1/4. It is seen that this depen- 
dence approaches a power law. The fitted exponent 
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Fig. 3. Figure a shows, in a log-log scale, the retreat distance 
L(t) and its power-law fit 1.73 1 0,50 for the case of m = 1/4. 
Figure b shows, in a log-log scale, the local maximum inter- 
face elevation, h max (t), and its power-law fit 0.97i 0,25 . 
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Fig. 4. Self-similarity at m = 1/4. Shown is the shape 
function h(xi,t), rescaled to the local maximum elevation 
h max , versus the coordinate xi, rescaled to h max , at times 
t = 8.35 X 10 3 (the black solid line), 1.07 X 10 5 (the red 
dashed line) and 1.16 X 10 6 (the blue dotted line). 

is 0.25, in excellent agreement with the theoretical 
value 1/4. Figure 4 demonstrates the presence of a 
self-similar region in the shape function for m = 1 /4. 
Also noticeable is a rapid (in time) convergence to 
the self-similar shape in the lobe region, and a slower 
convergence in the neck (the neck can be identified 
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with region II of our theory, as in Ref. [12]). 

The numerical results for to = 1/2 are presented 
in Figures 5, 6 and 7. Here too power laws for L(t) 
and R(t) are observed, and the fitted exponents 0.44 
and 0.28 are in good agreement with theoretical val- 
ues 3/7 ~ 0.43 and 2/7 ~ 0.29, respectively. The 
shape function again shows self-similarity, with a 
rapid (in time) convergence in the lobe region, and 
a much slower convergence in the neck. 

200-1 1 1 . 1 1 1 . 1 . 1- 
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Fig. 5. Snapshots of a part of the simulated system for 
m = 1/2 at times t = (the black solid line), 1.01 X 10 6 (the 
red dashed line) and 1.1 x 10 7 (the blue dotted line). No- 
tice the large difference between the horizontal and vertical 
scales. 

The numerical results for m — 5/4 are shown in 
Figures 8, 9 and 10. Here the self-similar regime is 
observed at very short times. The observed expo- 
nents of the power laws (0.29 for a and 0.35 for j3 
are in good agreement with theoretical values a — 
0.3 and f3 = 0.35. The shape function shows self- 
similarity in the lobe region, and a rapid deterio- 
ration of self-similarity in the neck as the time in- 
creases. 

Our numerical results for the dynamic exponents 
a and (3 at different to are summarized in Fig. 11. 
It can be seen that they follow theoretical curves 
predicted by Eq. (16). The numerically found coef- 
ficients a and b of the power laws are presented in 
Fig. 12. 

Figure 13 depicts, on a single plot, a set of rescaled 
shape functions for four different values of to. Three 
of them: for to = 1/4, 1/2 and 5/4, were computed 
in the present work, they are the same as shown in 
Figures 4, 7 and 10, respectively. The shape func- 
tion for to — is taken from Ref. [12]. Remarkably, 
all the shape functions coincide in the lobe region. 
That is, although the retreat distance and the local 
maximum elevation of the interface depend on to, 
the rescaled interface shape is independent of m, as 




Fig. 6. Figure a shows, in a log-log scale, the retreat distance 
L(t) and its power-law fit 1.43 t 0A4 for the case of m = 1/2. 
Figure b shows, in a log-log scale, the local maximum inter- 
face elevation, h max (t), and its power-law fit 1.24 1 - 28 . 
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Fig. 7. Self-similarity at m = 1/2. Shown is the shape 
function h(xi,t), rescaled to the local maximum elevation 
h max , versus the coordinate xi, rescaled to h max , at times 
t = 1.57 X 10 4 (the black solid line), 10 6 (the red dashed 
line) and 1.06 X 10 7 (the blue dotted line). 

long as to ^ 1. Figure 13 shows it very clearly in the 
lobe region. We believe, however, that the rescaled 
shape functions actually coincide in the neck too, 
and that the self-similar shape function computed 
for to = in Ref. [12] (see Fig. 13) is valid for any 
to 7^ 1. Unfortunately, it is hard to prove this con- 
jecture numerically. We observed that, at to < 1, 
non-self-similar corrections to the self-similar solu- 



7 



2.0x1 0" 6 - 
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00 1.0x10" 5 2.0x10 s 3.0x1 0" 5 

X 

Fig. 8. Snapshots of a part of the simulated system for 
m = 5/4 at times t = (the black solid line), 3.2 X 1CT 19 
(the red dashed line), 2.5 X 10 — 18 (the blue dotted line), and 
8.7 X 10 — 18 (the green dash-dotted line). Notice the large 
difference between the horizontal and vertical scales. 




1ff *1 10 ' ' 1Q- 19 10' 1 ' 



Fig. 9. Figure a shows, in a log-log scale, the retreat distance 
L(t) and its power-law fit 1.1 1 23 for the case of m = 5/4. 
Figure b shows, also in a log-log-scale, the local maximum 
elevation h max (t) (the circles) and its power-law fit 1.6t 0,35 . 

tion in the neck region decay very slowly with time 
when m is close to 1. As a result, one should go to 
prohibitively long times (and prohibitively large nu- 
merical domain sizes) in order to reach the "pure" 
similarity regime in the neck. In its turn, at m > 
1 corrections to the self-similar solution grow very 
fast with time when m is close to 1, so one has to 
go to prohibitively small times in order to observe 
the "pure" self-similarity in the neck. This is hard to 



Fig. 10. Self-similarity at m = 5/4. Shown is the shape 
function h(xi,t), rescaled to the local maximum elevation 
h max , versus the coordinate xi, rescaled to h max , at times 
t = 1.1 X 1CT 19 (the black solid line), 4.9 X 10~ 18 (the red 
dashed line), and 1.1 X 10 -16 (the blue dotted line). 

achieve in view of the presence of the numerical cut- 
off at m > 1 , which requires sufficiently long wait- 
ing times before its influence on the scaling results 
becomes small. 

0.6 -k- — ■ 1 ■ 1 ■ 1 ■ '- 

0.5- \. 
* 0.4- ^\ 

0.3- --"""""^"^a^^- 

0.2^1 , , , , , - 

0.0 0.5 1.0 1.5 

m 

Fig. 11. The exponents a and /3 of the power laws (8) versus 
the parameter m. The solid and dashed lines show theoretical 
predictions [see Eq. (16)] for a and /3, respectively. The 
corresponding numerical results are shown by the squares 
and circles, respectively. The numerical values of a and 
for m = 1/4, 1/2 and 5/4 are computed in this work, while 
those for m = and m = 1 are taken from Refs. [12] and 
[13], respectively. 

Is there any connection between the universal 
scaled shape shown in Fig. 13 and the rescaled 
interface shapes of the wedge (m = 1) which de- 
pend on the wedge angle? To address this question, 
we simulated the relaxation dynamics of several 
wedges with small angles. Figure 14 shows self- 
similar shape functions for three different values of 
the wedge angle: 10°, 5° and 2°. Shown on the same 
graph is the shape function for m = 0. Remarkably, 
all the shape functions coincide in the lobe region. 



8 



3 



2- 

■ 

0-U ■ , ■ , ■ , ■ , ■ r- 
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Fig. 12. The numerically found coefficients a (squares) and b 
(circles) of the power laws (8) versus the parameter m. The 
values for m = 1/4, 1/2 and 5/4 are computed in this work, 
while those for m = are taken from Ref. [12]. 
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Fig. 13. Self-similar shape functions at four different values 
of m: m = (the black solid line), 1/4 (the red dashed line), 
1/2 (the blue dotted line) and 5/4 (the green symbols). 

Furthermore, the numerics strongly suggests that, 
as the wedge angle tends to zero, the shape func- 
tion approaches that for m — everywhere. That 
is, the observed universal shape function for m^l 
coincides with the shape function obtained in the 
zero-angle limit for m = 1. 



1.5- 
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Fig. 14. Self-similar shape functions for m = 1 at three dif- 
ferent values of the wedge angle: 10° (the green dash-dotted 
line), 5° (the blue dotted line), and 2° (the red dashed line). 
The black solid line shows the shape function for m = 0. 



hibit self-similar intermediate asymptotics: a late- 
time asymptotics for < m < 1 and an early-time 
asymptotics for m > 1. Our theoretical predictions 
for the dynamic exponents of the retreat distance 
and the local maximum elevation of the interface 
versus time, at different m, are in excellent agree- 
ment with numerical simulations. We have found 
that, for m ^ 1, the rescaled interface shape is uni- 
versal, that is, it does not depend on m in the sim- 
ilarity region. Remarkably, the same rescaled inter- 
face shape also emerges in the case of m = 1 (where 
the self-similarity is exact and holds for the whole in- 
terface at all times t > 0), in the limit of zero wedge 
angle. 

Future work should provide a more complete the- 
ory that would explain the surprising findings pre- 
sented here. We hope that this work (see also Refs. 
[12,13]) will facilitate experimental and further the- 
oretical efforts aimed at a better understanding of 
the "simple" unforced Hele-Shaw flow. 



4. Summary 



We have investigated the dynamics of relaxation, 
by surface tension, of a family of curved interfaces, 
dividing an inviscid and viscous fluids in a Hele- 
Shaw cell, and characterizable by a single exponent 
m. A stripe m — 0, a wedge m = 1 and a generic 
cusp m = 2 that appears after a pinch-off event, rep- 
resent particular cases of these more general shapes. 
Combining simple analytic arguments with a ro- 
bust numerical method, we have found that, for any 
m^l, the relaxation dynamics of the interfaces ex- 
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